
**************** Welfare: Plots for Alternative Policy Designs ************

************* content *************
	** 1) load matrices and fix data structure
	** 2) scale variables for use in plot 
	** 3) construct 3D plot, wire option
			** outcome: y-variable 
	** 4) construct surface with contour lines 
			** outcome: z-variable 

***********************************

capture log close 
clear all
set more off 

cd "R:\WSV2\TBu_AKe\Spatial_NEW"

capture mkdir Graphs
global graphs "R:\WSV2\TBu_AKe\Spatial_NEW\Graphs"

log using welfare_3D_plot_complete, replace 


cd "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare"

capture ssc install 3dplot


	** 20x20 matrix over sigma and kappa for all outcomes 
use E // energy savings 

gen n = _n // merge would re-order obs 

forvalues s = 1(1)20 {
	rename c`s' e`s' // rename 
}

merge 1:1 _rowname using V // merge with V (utility loss)

drop _merge 

sort n

forvalues s = 1(1)20 {
	rename c`s' v`s' // rename 
}


merge 1:1 _rowname using W // merge with W (sales weights)

drop _merge 

sort n

forvalues s = 1(1)20 {
	rename c`s' w`s' // rename 
}

merge 1:1 _rowname using P // merge with P (restriction: yes/no)

drop _merge 

sort n

forvalues s = 1(1)20 {
	rename c`s' p`s' // rename 
}

	** fix data structure 
gen kappa = _n
gen sigma = _n

gen E_sum_R = . 
gen U_ratio = . 

gen R_ratio = . 
gen P_ratio = . 

fillin kappa sigma // data are placed as matrix, need long-form with 3 variables  

	** fill outcome variable over all combinations 
forvalues s = 1(1)20 {

sort kappa v`s'
bysort kappa: carryforward v`s', replace

sort kappa e`s'
bysort kappa: carryforward e`s', replace

sort kappa w`s'
bysort kappa: carryforward w`s', replace

sort kappa p`s'
bysort kappa: carryforward p`s', replace

}

	** place column values into single variable for long format 
forvalues k = 1(1)20 {
	forvalues s = 1(1)20 {
		
	replace E_sum_R = e`s' if sigma == `s' & kappa == `k' 
	replace U_ratio = v`s' if sigma == `s' & kappa == `k' 
	
	replace R_ratio = w`s' if sigma == `s' & kappa == `k' 
	replace P_ratio = p`s' if sigma == `s' & kappa == `k' 
		}
	}

	
gen k_index = kappa 
gen s_index = sigma 
replace kappa = 3+ kappa*0.1-0.1
replace sigma = sigma*0.01-0.01

sort kappa sigma 

save welfare_20x20, replace 

*use Welfare_3D_0802, replace


	*** implement 3D: scale outcomes 
sum U_ratio 
gen U_ratio2 = U_ratio*10 // roughly on 0 to 6 scale 

sum E_sum_R 

scalar sc =  3073453/5.915121 // scale factor 

gen E_sum_2 = E_sum_R/`=sc'
sum E_sum_2 

sum P_ratio R_ratio 

gen P_ratio2 = P_ratio*6 // max to 6. 

gen R_ratio2 = R_ratio*6 

sort sigma kappa 

		*******************
		**** 3D Plots *****
		*******************
cd "$graphs" 
		
		** private marginal utility loss 
xtset s_index k_index // treat sigma as panel, kappa as time. 
gen V_sum_R = U_ratio*E_sum_R
gen lV = L1.V_sum_R 
gen dV = ln(V_sum_R) - ln(lV) // at group_id = 2: percentage change in V_sum_R as sigma increases from zero. 
gen dV2 = V_sum_R - lV 

gen lE = L1.E_sum_R 
gen dE = ln(E_sum_R) - ln(lE) // at group_id = 2: percentage change in V_sum_R aa sigma increases from zero. 
gen dE2 = E_sum_R - lE

gen d_VE_k = dV2/dE2
sum d_VE_k

gen x = 5.91/1.13 // scale factor
gen d_VE_2 = d_VE_k*x // scaling to get cube proportions 
sum d_VE_2

		*** Private marginal utility loss
graph3d s_index d_VE_2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(240) blv format("%12.0fc") equi markeroptions(msize(1.2)) 

graph save 	 3D_VE_wire.gph, replace
graph export 3D_VE_wire.eps, as(eps) fontface("Baskerville Old Face") replace
graph export 3D_VE_wire.png, as(png) width(1500) replace			
			
		*** Utility ratio 
graph3d s_index U_ratio2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(240) blv format("%12.0fc") equi markeroptions(msize(1.2))  

graph save 	 3D_V_wire.gph, replace
graph export 3D_V_wire.eps, as(eps) fontface("Baskerville Old Face") replace
graph export 3D_V_wire.png, as(png) width(1500) replace	

		*** Proportion of restricted sales 
graph3d s_index P_ratio2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(232) blv format("%12.0fc") markeroptions(msize(1.2))  

graph save 	 3D_R_wire.gph, replace
graph export 3D_R_wire.eps, as(eps) fontface("Baskerville Old Face") replace
graph export 3D_R_wire.png, as(png) width(1500) replace	


graph3d s_index P_ratio2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(210)  blv format("%12.0fc") markeroptions(msize(1.2))  

graph save 	 3D_R_wire210.gph, replace
graph export 3D_R_wire210.eps, as(eps) fontface("Baskerville Old Face") replace
graph export 3D_R_wire210.png, as(png) width(1500) replace	


		*** Energy Savings  
graph3d s_index E_sum_2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(235) blv format("%12.0fc") equi markeroptions(msize(1.2))  

graph save 	 3D_E_wire.gph, replace
graph export 3D_E_wire.eps, as(eps) fontface("Baskerville Old Face") replace
graph export 3D_E_wire.png, as(png) width(1500) replace	

	
	***********************
	**** Contour Plots ****
	***********************
	
/* 
clear 
use welfare_100x100, replace // 100x100 simulation 
*/ 
	
		*** Private marginal utility loss  
twoway contour d_VE_k sigma kappa, yscale(reverse) color(hue%60) ecolor(red) colorlines levels(50) ztitle("") zlabel(0 0.09 "ABR (0.09)" 0.5 "0.5" 0.68 "Flat (0.68)" 1.12 "Max (1.12)", format(%9.2f) labsize(small)) xlabel(3 4.407 "4.407" 5, labsize(small)) ylabel(0 0.126 "0.126" 0.20, labsize(small)) ytitle("") xtitle("") name(ctour_VE, replace) graphregion(color(white))

graph save 	 Contour_VE.gph, replace
graph export Contour_VE.eps, as(eps) fontface("Baskerville Old Face") replace
graph export Contour_VE.png, as(png) width(1500) replace	


		*** Utility ratio 
twoway contour U_ratio sigma kappa, yscale(reverse) color(hue%75) ecolor(red) colorlines ccuts(0(0.05)0.6) ztitle("Loss per Unit of Energy Savings")  zlabel(0 0.133 "ABR (0.13)" 0.20 0.352 "Flat (0.35)" 0.40 .60, format(%9.2f) labsize(small)) xlabel(3 4.407 "4.407" 5, labsize(small)) ylabel(0 0.126 "0.126" 0.20, labsize(small)) ytitle("") xtitle("") name(ctour_V, replace) graphregion(color(white))

graph save 	 Contour_V.gph, replace
graph export Contour_V.eps, as(eps) fontface("Baskerville Old Face") replace
graph export Contour_V.png, as(png) width(1500) replace	



 	*** Energy Savings 
twoway contour E_sum_R sigma kappa  if E_sum_R, yscale(reverse) color(hue%60) ecolor(red) colorlines levels(50) ztitle("Total Energy Savings") zlabel(0 90000 "ABR (39 K)" 1000000 "1.0 M" 1801995 "Flat (1.8 M)" 3000000 "3.0 M", format(%9.2f) labsize(small)) xlabel(3 4.407 "4.407" 5, labsize(small)) ylabel(0 0.126 "0.126" 0.20, labsize(small)) ytitle("") xtitle("") name(ctour_E, replace) graphregion(color(white))

graph save 	 Contour_E.gph, replace
graph export Contour_E.eps, as(eps) fontface("Baskerville Old Face") replace
graph export Contour_E.png, as(png) width(1500) replace	


	*** Proportion of affected Sales 
twoway contour R_ratio sigma kappa, yscale(reverse) color(hue%60) ecolor(red) colorlines levels(50) ztitle("Fraction of Sales under Restriction") zlabel(0 0.21 "ABR (0.21)" 0.5  1, format(%9.2f) labsize(small)) xlabel(3 4.407 "4.407" 5, labsize(small)) ylabel(0 0.126 "0.126" 0.20, labsize(small)) ytitle("") xtitle("") name(ctour_R, replace) graphregion(color(white))
	
graph save 	 Contour_R.gph, replace
graph export Contour_R.eps, as(eps) fontface("Baskerville Old Face") replace
graph export Contour_R.png, as(png) width(1500) replace	


log close 
clear 
exit 


		*** Aggregate Utility Loss [outtake] ***
/* 
		*** totals utility loss (sum Delta V)
gen V_sum_R = E_sum_R*U_ratio 

scalar sc =  1817984/5.915121 // scale factor 

gen V_sum_2 = V_sum_R/`=sc'
sum V_sum_2 


graph3d s_index V_sum_2 k_index, cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(210) blv format("%12.0fc") markeroptions(msize(1.2)) 

graph save "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\3D_Vsum.gph", replace 
graph export "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\3D_Vsum.png", as(png) replace


	*** contour 
twoway contour V_sum_R sigma kappa, yscale(reverse) color(hue%60) ecolor(red) colorlines levels(50) ztitle("Sum Delta V") zlabel(minmax 5400 "ABR (5.4k)", format(%9.2f) labsize(small)) xlabel(3 4.407 "4.407" 5, labsize(small)) ylabel(0 0.126 "0.126" 0.20, labsize(small)) ytitle("") xtitle("") name(ctour_Vsum, replace) graphregion(color(white))

graph save "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\Contour_Vsum.gph", replace 
graph export "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\Contour_Vsum.png", as(png) replace


			*** plot kappa for given sigma *** 
			
gen V_sum_1 = V_sum_2 / 6 
gen U_ratio1 = U_ratio2 / 6 
gen E_sum_1 = E_sum_2 / 6 


			
line V_sum_1 U_ratio1 E_sum_1 kappa if s_index == 1, xline(4.407) name(comp_EUV00, replace) title("sigma = 0") yscale(range(0 1))

line V_sum_1 U_ratio1 E_sum_1 kappa if s_index == 8, xline(4.407) name(comp_EUV08, replace) title("sigma = 0.08") yscale(range(0 1))

line V_sum_1 U_ratio1 E_sum_1 kappa if s_index == 12, xline(4.407) name(comp_EUV12, replace) title("sigma = 0.12") yscale(range(0 1))

line V_sum_1 U_ratio1 E_sum_1 kappa if s_index == 16, xline(4.407) name(comp_EUV16, replace) title("sigma = 0.16") yscale(range(0 1))


grc1leg comp_EUV00 comp_EUV08 comp_EUV12 comp_EUV16, name(comp_EUVall, replace) legendfrom(comp_EUV00)

graph save "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\comp_EUVall.gph", replace 
graph export "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\comp_EUVall.png", as(png) replace


line V_sum_1 U_ratio1 E_sum_1 sigma  if k_index == 1, xline(0.126) name(comp_EUV3, replace) title("kappa = 3") yscale(range(0 1)) xscale(reverse)

line V_sum_1 U_ratio1 E_sum_1 sigma  if k_index == 10, xline(0.126) name(comp_EUV4, replace) title("kappa = 4") yscale(range(0 1)) xscale(reverse)

line V_sum_1 U_ratio1 E_sum_1 sigma  if k_index == 15, xline(0.126) name(comp_EUV44, replace) title("kappa = 4.4") yscale(range(0 1)) xscale(reverse)

line V_sum_1 U_ratio1 E_sum_1 sigma  if k_index == 20, xline(0.126) name(comp_EUV5, replace) title("kappa = 5") yscale(range(0 1)) xscale(reverse)


grc1leg comp_EUV3 comp_EUV4 comp_EUV44 comp_EUV5, name(comp_EUVsig, replace) legendfrom(comp_EUV3)

graph save "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\comp_EUVsig.gph", replace 
graph export "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\comp_EUVsig.png", as(png) replace


scalar sigma = 0 
scalar kappa = 4.407 
scalar gamma = 0.126



twoway (connected V_sum_R E_sum_R if s_index == 1) (connected V_sum_R E_sum_R if s_index == 4) (connected V_sum_R E_sum_R if s_index == 8) (connected V_sum_R E_sum_R if s_index == 12) (connected V_sum_R E_sum_R if s_index == 16)

twoway (connected V_sum_2 E_sum_2 if s_index == 1) (connected V_sum_2 E_sum_2 if s_index == 4) (connected V_sum_2 E_sum_2 if s_index == 8) (connected V_sum_2 E_sum_2 if s_index == 12) (connected V_sum_2 E_sum_2 if s_index == 16)

replace E_sum_1 = E_sum_1*20

graph3d s_index V_sum_2 E_sum_1 , cuboid innergrid wire colorscheme(bcgyr) xlab("") ylab("") zlab("") xlang(10) ylang(90) zlang(310) xlpos(2) ylpos(2) zlpos(9) xang(12) yang(250) blv format("%12.0fc") markeroptions(msize(1.2)) 

graph save "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\3D_Vsum.gph", replace 
graph export "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare\3D_Vsum.png", as(png) replace

*/ 

